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Abstract. Photoacoustic computed tomography (PACT), also known as optoacous- 
tic tomography, is an emerging imaging modality that has great potential for a wide 
range of biomedical imaging applications. In this Note, we derive a hybrid recon- 
struction formula that is mathematically exact and operates on a data function that 
is expressed in the temporal frequency and spatial domains. This formula explicitly 
reveals new insights into how the spatial frequency components of the sought-after 
object function are determined by the temporal frequency components of the data 
function measured with a circular or spherical measurement geometry in two- and 
three-dimensional implementations of PACT, respectively. The structure of the re- 
construction formula is surprisingly simple compared with existing Fourier-domain re- 
construction formulae. It also yields a straightforward numerical implementation that 
is robust and two orders of magnitude more computationally efficient than filtered 
backprojection algorithms. 
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Photoacoustic computed tomography (PACT) is an emerging imaging modality 
that has great potential for a wide range of biomedical imaging applications (Oraevsky 
& Karabutov, 2003; Wang, 2008; Kruger, et al., 1999). In PACT, biological tissues of 
interest are illuminated by use of short laser pulses, which results in the generation 
of internal acoustic wavefields via the thermoacoustic effect (Xu & Wang, 2006; Xu, 
et al., 2010). The initial amplitudes of the induced acoustic wavefields are proportional 
to the spatially variant absorbed optical energy density within the tissues. The 
propagated acoustic wavefields are subsequently detected by use of a collection of 
ultrasonic transducers that are located outside the object. An image reconstruction 
algorithm is employed to estimate the absorbed optical energy density within the tissue 
from these data. 

A variety of image reconstruction algorithms have been proposed for PACT (Xu, 
et al, 2002; Finch, et al, 2004; Xu & Wang, 2005; Finch, et al, 2007; Kunyansky, 
2007; Hristova, et al., 2008; Kunyansky, 2012; Salehin & Abhayapala, 2012; Treeby 
& Cox, 2010). While iterative image reconstruction methods hold great value due to 
their ability to incorporate accurate models of the imaging physics and the instrument 
response (Paltauf, et al., 2002; Yuan & Jiang, 2007; Zhang, et al., 2009; Provost 
& Lesage, 2009; Wang, et al, 2011; Guo, et al., 2010; Huang, et al, 2010; Xu, 
et al, 2011; Buehler, et al, 2011; Bu, et al., 2012; Wang, et al, 2012), they can lead to 
long reconstruction times, even when accelerated by use of modern computing hardware 
such as graphics processing units (Wang et al., 2012). This is especially problematic in 
three-dimensional (3D) implementations of PACT, in which reconstruction times can 
be excessively long. Almost all experimental studies of PACT to date have employed 
analytic image reconstruction algorithms. Even if an iterative image reconstruction 
algorithm is to be employed, it is often useful to employ an analytic reconstruction 
algorithm to obtain a preliminary image that can initialize the iterative algorithm and 
thereby accelerate its convergence. 

Most analytic reconstruction algorithms for PACT with a spherical measurement 
aperture and point-like transducers have been formulated in the form of filtered 
backprojection (FBP) algorithms. These algorithms possess a large computational 
burden, requiring 0(N 5 ) floating point operations to reconstruct a 3D image of 
dimension iV 3 . Image reconstruction algorithms based on the time-reversal principle 
and finite-difference schemes require 0(iV 4 ) operations (Burgholzer, et al., 2007). 
Fast reconstruction algorithms for spherical measurement apertures that require only 
0(N 3 logN) operations have been proposed (Kunyansky, 2012; Salehin & Abhayapala, 
2012). However, numerical implementations of these formulas require computation 
of special functions and multidimensional interpolation operations in Fourier space, 
which require special care to avoid degradation in reconstructed image accuracy. It 
is well-known that the temporal frequency components of the pressure data recorded 
on a spherical surface are related to the Fourier components of the sought-after object 
function (Anastasio, et al., 2007). However, to date, a simple reconstruction algorithm 
based on this relationship, i.e., one that does not require series expansions involving 
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special functions or multi-dimensional interpolations, has yet to be developed. 

In this Note, we derive a novel reconstruction formula for two-dimensional (2D) and 
3D PACT employing circular and spherical measurement geometries, respectively. The 
mathematical forms of the reconstruction formulae are the same in both dimensions and 
are surprisingly simple compared with existing Fourier-domain reconstruction formulae 
for spherical and circular measurement geometries. The reconstruction formulae are 
mathematically exact and describe explicitly how the spatial frequency components of 
the sought-after object function are determined by the temporal frequency components 
of the measured pressure data. Their discrete implementations require only discrete 
Fourier transform, one-dimensional interpolation, and summation operations. A 
preliminary computer-simulation study is conducted to corroborate the validity of the 
reconstruction formula. 

We consider the canonical PACT imaging model in which the object and 
surrounding medium are assumed to possess homogeneous and lossless acoustic 
properties and the object is illuminated by a laser pulse with negligible temporal 
width. Point-like, unfocused, ultrasonic transducers are assumed. We also assume 
that the effects of the acousto-electric impulse responses of the transducers have been 
deconvolved from the measured voltage signals so that the measured data can be 
interpreted as pressure signals. The 3D problem is addressed where p(r, t) denotes 
the photoacoustically-induced pressure wavefield at location r G M 3 and time t > 0. 
However, the analysis and reconstruction formula that follows remains valid for the 2D 
case. The imaging physics is described by the photoacoustic wave equation (Oraevsky 
& Karabutov, 2003; Wang, 2008; Kruger et al, 1999): 

„9 / 1 d 2 p(r,t) . . 

subject to the initial conditions: 

P (r, f ) =^(r); =0, (2) 

*=o Up at t=o 

where V 2 denotes the 3D Laplacian operator and A(r) is the object function to be 
reconstructed that is contained within the volume V. Physically, A(r) represents the 
distribution of absorbed optical energy density. The constant quantities j3, c, and C p 
denote the thermal coefficient of volume expansion, speed-of-sound, and the specific 
heat capacity of the medium at constant pressure, respectively. 

Let p(r s ,t) denote the pressure data recorded at location r s G S on a spherical 
surface S of radius Rs that encloses V. The continuous form of the imaging model 
that relates the measurement data to object function can be expressed as (Cox & 
Beard, 2005): 

p(r„ t)=UA = [ dki(k) cos(cH)e Jk - r % (3) 

U p {ZTT) Joe 

where k 6 I 3 is the spatial frequency vector conjugate to r, k = |k|, and A(k) is the 
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3D Fourier transform of A(r). We adopt the Fourier transform convention 

i(k) = T 3 A(r) = [ drA(v)e-^ r (4a) 



A(r) = ^A(k) ee JL £ dki(k)e ?kr . (46) 

The imaging model in Eqn. fl3]) can be interpreted as a mapping "H : O — > D between 
infinite dimensional vector spaces that contain the object and data functions. We will 
define O as the vector space of bounded and smooth functions that are compactly 
supported within the volume V . 

Let the infinite set of functions {7 M (r)}, indexed by fi, represent an orthonormal 
basis for O. The object function A(r) can be represented as 



A(r)= / dniA^Mr), (5) 

J oo 

where the inner product in O is defined as 

(A, 7 „) = £drA(r) 7 ,(r) = 7^/ dkA(k)%(k), (6) 

7 M (k) = J r 3 7 M (r), and the quantity on the right-hand side of Eqn. ([SD follows from fact 
that the Fourier transform is an isometry. A trace identity (see Eqn. (1.7) in reference 
(Finch et al, 2004) for the 3D case and Eqn. (1.16) in (Finch et al., 2007) for the 2D 
case) can be employed to relate the inner products in the spaces O and D as: 

2C^ f°° f 

(^^) = 1T-^/ dt I dx 8 tp{T a ,t)v^T 8 ,t), (7) 
tt-SP Ct Jo Js 

where 

v„(r„ t) = = J dk 7 (k) cos(cH)e Jk - (8) 

and the right-hand side of Eqn. (jj]) defines a scaled version of the inner product in D. 
On substitution from Eqn. (jSJ) into Eqn. ([7j), one obtains 

(A7m) = 7^)3 / ^k^(k) 7M (k), (9) 

where 

y(k) = — ± / dr s e ik ^ / dttp(r s , t) cos(ctt). (10) 
ttsP Js Jo 

Comparison of Eqns. ([6]) and (Q reveals that A(k) = y(k). By evaluating the Fourier 
cosine transform that is present in the right-hand side of Eqn. ( TlOl) . a reconstruction 
formula for determining A(k) can therefore be expressed as 

i(k) = ^/^ e ' k ' rsRe { J ' l{ ^ (r '- t)}(rs ' a;)L =^}' (U) 
where T\ denotes the one-dimensional (ID) Fourier transform with respect to time t 
and 'Re' denotes the operation that takes the real part of quantity in the brackets. 
Subsequently, A(r) is determined as J r 3 " 1 A(k). 
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Equation (II ip represents a novel reconstruction for PACT and is the key result 
of this Note. Unlike previously proposed Fourier- domain reconstruction formulae 
(Norton, 1980; Kunyansky, 2012; Salehin & Abhayapala, 2012), Eqn. §Tl§ has a 
simple form and does not involve series expansions utilizing special functions. The 
reconstruction formula reveals that the measured data p(r s , t) determine the 3D Fourier 
components of the A(r) via a simple process that involves the following four steps: (1) 
Compute the ID temporal Fourier transform of the modified data function tp(r s , £); (2) 
Isolate the real- valued component of this quantity corresponding to temporal frequency 
u = ck; (3) Weight this value by the plane-wave e ikrs ; and (4) Sum the contributions, 
formed in this way, corresponding to every measurement location r s £ S. This reveals 
the components of A(k) residing on a sphere of radius - are determined by the ID 
Fourier transform of tp(r s , t) corresponding to temporal frequency u. In this sense, Eqn. 
( ITT]) can be interpreted as an implementation of the Fourier Shell Identity (Anastasio 
et al., 2007). Finally, the form of Eqn. (TTTj) remains unchanged in the 2D case, where 
r s , k £ M 2 and S is a circle that encloses the object. 

A discrete implementation of Eqn. (ITT)) possesses low computational complexity and 
desirable numerical properties. The ID fast Fourier transform (FFT) can be employed 
to approximate the action of T\ and only a ID interpolation is required to determine 
the value of the Fourier transformed data function corresponding to temporal frequency 
uj = ck, where k corresponds to the magnitude of vectors k that specify a 3D Cartesian 
grid. From the values of A(k) determined on this grid, the 3D FFT algorithm can be 
employed to estimate values of A(r). If the object is represented on a N x N x N 
grid and the number of transducer locations and time samples are both O(N), the 
computational complexity is limited by the 3D FFT algorithm, i.e., 0(N 2 \ogN) in 2D 
and O(iV 3 logA0 in 3D. 

A preliminary computer-simulation study for the 2D case was conducted to 
corroborate the correctness of the reconstruction formula. The object function A(r) 
was taken to be the numerical phantom shown in Fig. [H-(a), which was comprised of 
a collection of uniform disks that were blurred by a 2D Gaussian kernel whose full- 
width-at-half-maximum was 0.3 mm. The phantom was discretized by use of pixel 
expansion functions with a pitch of 0.025 mm. The measurement geometry consisted 
of 256 point-like transducers that were uniformly distributed over a circle of radius 
12.8 mm that enclosed the object. The k-wave toolbox (Treeby & Cox, 2010) was 
employed to numerically solve the photoacoustic wave equation and generate simulated 
pressure signals at each transducer location at a temporal sampling rate of 30 MHz. The 
simulated pressure data set generated in this way contained 256 x 2048 data samples. 
The speed-of-sound and -jj- were assigned values of 1.5-mm//is and 1000 (arbitrary 
units), respectively. A noisy data set was produced by addition of 5% uncorrelated 
Gaussian noise to the noiseless pressure data. 

Images were reconstructed on a uniform 2D grid of spacing 0.1 mm by use of a 
discretized form of Eqn. (jll|) coupled with the 2D inverse FFT algorithm. In order 
to reconstruct images of dimension 256 x 256, samples of A(k) were determined on a 
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uniform 2D grid of dimension 512 x 512 with a sampling interval of (0.1 x 256) -1 mm -1 . 
The samples of the data function tp(r s ,t) were zero-padded by a factor of 8 prior to 
estimating its ID Fourier transform by use of the FFT algorithm. From these data, 
nearest neighbor ID interpolation was employed to determine the values of the term in 
brackets in Eqn. fTTTT) corresponding to u = ck for the sampled locations k. 

The images reconstructed from the noiseless and noisy data sets are shown in Fig. [I]- 
(b) and (c). Profiles corresponding to the central rows of these images are shown in Fig. 
[2j These results confirm that the proposed reconstruction algorithm can reconstruct 
images with high fidelity from noise-free measurement data. Although, a systematic 
investigation of the noise propagation properties of the proposed algorithm is beyond 
the scope of this Note, Figs. [Er(c) and [2]-(b) suggest that its performance is robust 
in the presence of noise. This is to be expected, since all operations involved in the 
implementation of Eqn. ( ITT]) are numerically stable. 

In summary, we have derived a Fourier-based reconstruction formula for PACT 
employing circular and spherical measurement apertures. The formula is mathematically 
exact and possesses a surprisingly simple form compared with existing Fourier- 
domain reconstruction formulae. The formula yields a straightforward numerical 
implementation that is stable and is two orders of magnitude more computationally 
efficient than 3D filtered backprojection algorithms. The proposed formula serves as 
an alternative to existing fast Fourier-based reconstruction formulae. A systematic 
comparison of the proposed reconstruction formula with existing formulae by use of 
experimental data remains an important topic for future studies. 
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Figure 1. The numerical phantom is shown in subfigure (a). Images reconstructed by 
use of the proposed reconstruction algorithm from noiseless and noisy data are shown 
in subfigures (b) and (c), respectively. The greyscale window is [—0.2, 1.2]. 
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Figure 2. Profiles corresponding to the central rows of the images shown in Fig. Q} 
(b) (subfigure(a)) and Fig. [T]-(c) (subfigure(b)). The solid line in subfigure (a), which 
corresponds to the image reconstructed from noiseless data, almost completely overlaps 
with the profile through the numerical phantom. 



